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Shear flow controlled mode selection in a nonlinear autocatalytic medium 
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The effect of shear flow on mode selection and the length scale of patterns formed in a nonlinear auto-catalytic 
reaction-diffusion model is investigated. We predict analytically the existence of transverse and longitudinal 
modes. The type of the selected mode strongly depends on the difference in the flow rates of the participating 
species, quantified by the differential flow parameter. New spatial structures are obtained by varying the length 
scale of individual modes and superposing them via the differential flow parameter. Our predictions are in line 
with numerical results obtained from lattice Boltzmann simulations. 
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I. INTRODUCTION 


Reaction-diffusion models of excitable media are known to 
give rise to spatially and/or temporally varying patterns under 
certain conditions. The formation of these patterns in station¬ 
ary media can be explained based on a mechanism first pro¬ 
posed by Turing iQl]. However, in moving excitable media, the 
mechanism of pattern formation, length scale of the patterns 
formed and the type of patterns selected can be significantly 
different from the stationary case. Under a constant uniform 
flow for instance, differences in the flow rates of the partic¬ 
ipating reacting species can spatially disengage the species 
leading to the well known differential flow instability (DIFI) 
mechanisms Min. Patterns formed from this mechanism 
are known to be entirely different from those resulting from 
the Turing mechanism occurring in stationary media. In a lin¬ 
ear shear flow, the effect of the distribution of fluid velocities 
has been shown to give rise to another mechanism of pattern 
formation which does not require differences in flow rate of 
the reacting species or the fulfillment of Turing condition |@- 
1^. The assumption of a reacting fluid comoving with the av¬ 
erage flow velocity is not sufficient to explain these effects. A 
variety of different approaches have been proposed to study 
this problem. In Ref. lUOr . for example, a system of two cou¬ 
pled layers, moving at a constant velocity with respect to one 
another, is considered. This approach reproduces some of the 
key features of pattern formation in a shear flow including sit¬ 
uations where the reactants move at different flow rates but 
with the same diffusion rate d. However, these studies are 
performed for a piecewise defined velocity and differential 
flow rates for a ID model and do not include the combined 
effects of shear, differential advection and differential diffu- 
sivity of the reacting species on pattern selection. The roles 
played by both differential advection and diffusion of the re¬ 
acting species on pattern selection can often not be separated. 
As we have recently shown EQ, such a situation may be 
relevant when modelin g, e .g., growth kinetics and spatial dis¬ 
tribution of vegetation 1114 In this case, one of the react- 
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ing species (water) undergoes convection while the other com¬ 
ponent (vegetation) is transported via diffusion only. Other 
examples here are interaction of a flowing chemical species 
with another chemical species bound to a catalytic surface as 
in packed bed reactors lUa [O or cell polarization in devel¬ 
opmental biology S. 

In this work, we address this issue both analytically as well as 
via 2D numerical simulation of the nonlinear equations. Via 
a linear stability analysis of the advection-diffusion-reaction 
equations, we first determine the dispersion relation for the 
growth rate of perturbations for an arbitrary differential flow 
parameter. A simple scaling ansatz then allows to also include 
the effect of shear rate in the characteristic length scale of the 
patterns. We And that, by tuning the differential flow parame¬ 
ter, it is possible to select between transverse, longitudinal or 
mixed modes. In the latter case, the interaction of the selected 
modes can give rise to interesting novel structures, where the 
underlying length scale is tuned by shear rate. 


II. THE MODEL 


As a prototypical example, we choose the two species Gray- 
Scott model with an advection term to illustrate the effect of 
shear flow on pattern formation. The most general situation, 
as observed in aqueous solutions of reactants, is the coupling 
of the nonlinear reaction to the ambient fluid flow via the den¬ 
sity or viscosity of the fluid |T^ Such situations lead to 
hydrodynamic instabilities and fluxes in the horizontal layer 
of the reacting fluid ll2ll] . Rather than considering this case, 
we focus here on the simple situation of the so called ’passive 
advection’, where the flow field is imposed externally and de¬ 
coupled from the reacting species. The resulting advection- 
diffusion-reaction (ADR) equation describing the transport of 
the two species A and B can then be written as ifl^ 
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(l-A)-B^A + rADAV^A, ( 1 ) 
—B A- rjB^A + tbDb'V^B, (2) 


where A = A{x,y,t) and B = B{x,y,t) are the dimen¬ 
sionless concentrations of the interacting species. The time 
scales TA and tb are the characteristic time for the addition 
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and removal of the species A and B respectively. They are 
related to the rate constants kf and k 2 as, = 1/kf and 
tb = l/{kf + k2) lO- The parameter p sets the activation of 
the interacting species and is related to the reaction rate con¬ 
stants fci as p = Ao{kikf)^^'^/{kf + ^ 2 ) ifl^ . Da and Db 
are the diffusion coefficient of the two species and u{x, y) is 
the local flow velocity. The parameter 5, designated here as 
the differential flow parameter, is the ratio of the flow rates of 
the two interacting species. For simplicity, we consider here 
the flow between two parallel walls driven either by gravity 
(Poiseuille) or by the relative motion of the walls with respect 
to one another (Couette). For the Poiseuille flow, the velocity 
profile reads u{y) = u{y)x = {du^vg/Ly)y{Ly — y)x. Here, 
X is the unit vector along the x direction, Uavg is the average 
fluid velocity and Ly denotes the channel width. For the pla¬ 
nar Couette flow, u{y) = '^yx, with a constant shear rate, 
7- 

In the absence of advection, the system in Eqs. ([TJ and B ad - 
mits three spatially homogeneous steady state solutions 111 ^ . 
The first solution is the trivial homogeneous state Bf. = 0, 
Ae = 1 and exists for all system parameters. The other two 
solutions exist provided that the parameter 77 > 2. They are 
given by 
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and 
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(3) 

These solutions are destabilized by Turing and/or Hopf in¬ 
stability, leading to spatially and/or temporally varying struc¬ 
tures. The conditions for Turing and Hopf instabilities in the 
system are Da/ Dr > {tb/TA) r]B~ and ta/tb > r]B~, 
respectively fo . 


III. RESULTS 
A. Numerical simulation 

A first glance of flow effects on the resulting patterns is ob¬ 
tained by solving Eqs. ([U and B numerically in a rectangular 
domain with an imposed Poiseuille flow along the a;-direction. 
The numerical solution is performed with the lattice Boltz¬ 
mann method ll^ . Via a multiscale expansion analysis, we 
have shown in a previous work that, adding appropriate re¬ 
action terms to the standard lattice Boltzmann method allows 
to recover Eqs. O and B ifT^ . The channel dimension is 
100 X 30 lattice units. Periodic boundary condition is ap¬ 
plied along the flow direction. At the walls of the channel, 
a no-flux boundary condition is imposed for the concentration 
field, while the no-slip condition is applied to the fluid veloc¬ 
ity. The initial condition consists of a random perturbation 
to the steady state solution Eq. B at the channel inlet. Note 
that this is different from a fixed profile as used in The 
diffusion coefficients of the two species are chosen such that 
the ratio Da/Db satisfies the condition for Turing instabil¬ 
ity iQ. 

Results obtained from these simulations are depicted in Eig. 
[T] We observe longitudinal and transverse modes and super¬ 



FIG. 1: (Color online) Patterns obtained from numerical simulations 
of the model at different differential flow parameters, (a) 5 = 0, (b) 
5 = 1, (c) 5 = 0.7 [excitation of both the modes in (a) and (b)]. (d) 
5 = 0.7[consisting of the mode in (a) and a mode similar to (b) but 
containing a single stripe]. The system parameters are chosen as 
77 = 2.1988,I7A/f/B = '&.0,ta/tb = 2.744 and 7 = 0.001 and 
are the same in each panel except in (d) where Da/Db ~ 4.0. 


position thereof depending on the value of the differential flow 
parameter. Note that these patterns are stable and steady over 
time compared to the time for the slowly diffusing species to 
sample the width of the channel (i.e. for t ~ L^/Db)- Here, 
“transverse” and “longitudinal” refer to the orientation of the 
wavefront with respect to the flow direction. If the wavefront 
is perpendicular (parallel) to the flow direction then the mode 
is called transverse (longitudinal). Above a critical shear rate, 
advecting one of the species, while keeping the other immo¬ 
bile (S = 0) leads to a transverse mode with wave vector 
Q = {Qx, 0) (Eig. [Ha)). Eurther increase in the shear rate does 
not change the type of the mode selected, it however leads to 
the change in wavelength of the pattern. Advecting the two 
species at the same flow rate (6 = 1 ) leads to longitudinal 
modes with wave vector q ~ ( 0 , qy) (Eig.Bb))- In this case, 
both the type of the selected mode and its wavelength are in¬ 
dependent of shear rate. Eor 1 < 5 < 0, both the longitudinal 
and transverse modes are exited, giving rise to interesting spa¬ 
tial structures as exemplified in Eigs.[TJc) and (d). EigureBc); 
for instance, is obtained from the interaction of the four stripe 
transverse mode in (a) with the two stripe longitudinal mode in 
(b), while the structure in (d) is obtained from the interaction 
of the four stripe transverse mode in (a) with a single-stripe 
longitudinal mode. 


B. Shear dispersion effects 

The structures observed above in Eig. [1] are due to the in¬ 
homogeneity of the flow foil . To rationalize the dispersion 
effect introduced by the inhomogeneous flow field, we note 
that the characteristic time scales for the addition and removal 
of the species A and B are ta and tb- Within these time 
scales, species diffuse over length scales Ia = {2DataY^'^ 
and Ib = {^DbtbY^"^■ We have shown in previous work 
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that these length scales determine the characteristic length of 
the resulting patterns in the absence of flow d. Here we 
go one step further and show that modification of these length 
scales by the flow provides a simple way of tuning the char¬ 
acteristic length of the pattern along the flow direction and 
incorporating the effect of shear flow. For this purpose, we 
consider thermal diffusion of a particle of species A across 
streamlines of different velocities. Within a time scale of the 
order of ta, the particle diffuses a distance of the order of 
Ayo = V‘^Da Ta and experiences a change in the flow ve¬ 
locity of the order of Am = jAyQ. Its displacement along 
the direction of flow thus contains both a contribution due to 
thermal motion, Axq, and from the change of flow velocity, 
Ax^a” = taAua = ta'tAj/o- One can, therefore, write 

{if^Af = ((Axo + Ax«™)^) 

= 2DATA + iTAi?2DATA = 2Df,^TA (4) 

^'f,A = + {taiY ) ( 5 ) 

where we used the isotropy of thermal diffusion and the fact 
that the displacements due to the flow and diffusion are not 
correlated (AxoAcc*'”'"') = 0. The last equality in Eq. (|4]i de¬ 
fines the flow-enhanced effective diffusion coefficient along 
the flow direction = Da [l + ■ Similarly, and 

noting that Aub = jAyoS, one obtains for the species B, 


lf,B = lBV^ + irBiSy (6) 

and = Db [l -f (tb 7 i 5 )^] . It is noteworthy that our re¬ 
sults differ from the well-known Taylor dispersion by the fact 
that effective diffusion in the former continuously increases 
with time = D[1 + ( 7 ^)^]) 13, while in our 

case, characteristic time scales ta and tb set a limit to this 
process. 

The validity of Eq. (|5]l is examined by choosing b = 0, and 
performing series of simulations for two different flow situ¬ 
ations; a planar Couette flow and a Poiseuille flow between 
two planar plates. The effective characteristic length scale 
IA is obtained from the simulation results by computing the 
characteristic wavelength A for each pattern. As shown by 
linear analysis of the ADR equations for a constant uniform 
flow lfisll . the wavelength of the patterns is related to the char¬ 
acteristic length scale I a via {q = 2tt/X) 


^ ntA? 

Investigating the case (5 = 0 allows one to focus on the mod¬ 
ulation of species A, while species B is immobilized = 
Ib)- As seen from Eig.|2] the simple scaling theory accounts 
quite well for the variation of the characteristic length of the 
pattern with the flow. It is noteworthy that, while shear rate 
is constant across the channel in the first case, it varies lin¬ 
early with distance from the channel center in the case of the 
Poiseuille flow. Here, the shear rate appearing in Eq. Q is 
estimated by taking the ratio of the average flow velocity to 



FIG. 2: (Color online) Simulation results on the characteristic pat¬ 
tern length along the flow direction, (symbols), versus the pre¬ 
diction of the simple scaling theory, Eq. (S (solid lines). The simu¬ 
lation parameters are chosen as (5 = 0, Ia/Ib = 4.75, ta = 3030, 
Tb = 1104 and y = 2.1988. The initial configuration in all cases 
shown in the panel correspond to random perturbations added to the 
homogeneous steady state Eq. (O at the inlet of the channel. Panel 

(a) shows the result obtained for a simple shear flow (planar Couette 
geometry), where the shear rate is constant across the channel. Panel 

(b) corresponds to a parabolic velocity profile between two parallel 
plates (Poiseuille flow). Note that, in this case, the shear rate varies 
across the channel. Therefore, we use the average shear rate (defined 
as the mean flow velocity divided by the plate separation) when eval¬ 
uating Eq. Obviously, this simple scaling result provides a good 
approximation to the simulated data. 


the plate separation, 7 = u^vg/H. The good agreement be¬ 
tween theory and simulation in this case (Eig. |3) underlines 
the robustness of the scaling analysis, leading to Eq. ©. 
While tuning the magnitude of flow allows to set the length 
scale of the patterns, as already shown in Eig. [T] differential 
flow parameter, S, provides a mean to control mode selection. 
In the light of the above results, this observation can be put 
on a more rational basis. At first order, the flow modulates 
the diffusion coefficient in the flow (x-) direction while in the 
y-direction it remains unchanged, Eq. (HJi. Incorporating this 
effect and non-dimensionalizing time and length in Eqs. O 
and (| 2 ]i with a characteristic time scale ta and length scale I a, 
one can rewrite these equations as 


dA 

'm 


dA 

dx 


IdB 5 dB 
T dt T dx 


(1-41)- 

Hi^Af 


- B^A + V^A 
d'^A 


-B + qB^A+^V^B 

[S'^tbY d’^B 
£2 dx^ ’ 


( 8 ) 

(9) 


where r = ta/tb and e = Ia/Ib- 


C. Linear Stability Analysis 

In order to determine the selected mode, we perform a lin¬ 
ear stability analysis around the non-trivial homogeneous state 
Eq. (O and set 




















( 10 ) 

Here, (j)A -C Ag and (j)B -C Se are perturbation amplitudes 
and a is the growth rate of the perturbations. Note that the 
wavevectors are in dimensionless units such that = QxIa^ 
we have dropped the tilde symbol in subsequent discussion for 
clarity. Inserting this ansatz into Eqs. I© and ©, linearizing 
and requiring a non-trivial solutions for a leads to a charac¬ 
teristic equation for the growth rate of type 

+ {a + ib)a -f (c -I- id) = 0, (11) 

The coefficients a, b, c and d in Eq. (fTTIi are functions of the 
system parameters and the wave vector components Qx and qy. 
The coefficients are given as 


a = -q^{l + T/e^) + {T-r]Bt) 

-QliiT-A)^{S/e^ + 1), 

b = QxUavg{^5 1), 

c = rq^je^ ^ q^{rr]B^je‘^ - T - Sul^^) 

- 2) + Tq^{jTA)^/e^ 
+'r<lli'iTA)'^iSriBf/e'^ - 1) 
+{Tqlq'^{iTA)'^/£'^){S + 1 ), 
d = u^^g{q^^'^{T/e'^ + 6) + q{6riB^ - t) 

+qqU'iTA)^iS + St)), ( 12 ) 

where (g^ = ql + q^)- 

Evaluating the roots of the polynomial in Eq. (fTTI) and separat¬ 
ing the solution into the real and imaginary parts one obtains 


Re [a] 
Im[a] 


a ^ 1 Ir + (a^ — b^ — 4c) 

2 ^ 2 V 2 

6^1 Ir — [a? — b^ — 4c) 

2 ^ 2 V 2 


(13) 

(14) 


where r = (a^ — b^ — 4c)^ -f {2ab — 

A non-zero imaginary part of a indicates the presence of os¬ 
cillations, while the real part of a tells us whether an infinitesi¬ 
mal perturbation will decay (Re [a] < 0) or grow (Re [a] > 0), 
the latter being of particular interest for pattern formation. 
Moreover, the mode that has the largest positive real part of 
a is most probable to be selected. 

Figure |3^ shows the effect of shear rate on Re [a] as obtained 
from an evaluation of Eq. (fTTI) . In line with our scaling analy¬ 
sis of shear-induced diffusion, the wave vector associated with 
the fastest growing mode (maximum of Re [a]) decreases with 
increasing shear rate. Another observation from Fig.is that 
for a given value of e, the system becomes unstable upon in¬ 
creasing shear rate (see 7 > 0 .00015). An imposed shear flow 
thus has two effects; (i) it gives rise to new instabilities and 
(ii) it influences the length scale of the fastest growing mode. 
Next we turn to the effect of the differential flow parame¬ 
ter, 6, on the selected mode. For this purpose, we recall that 



FIG. 3: (Color online) (a) Effect of shear rate on the growth rate 
of the qx mode. The maximum of Re [a] shear rate showing a shift 
in the length scale of the patterns towards longer wavelength. The 
parameters are chosen as = 0, 5 = 0, e = 4.75, rj = 2.1988. 
(b) Effect of e on the growth rate of the qy mode. For e < rjBs 
[rjBf, = 3.421), the qy mode is damped. The parameters are chosen 
as = 0, (5 = 1.0, e = 4.75, p = 2.1988. 


diffusion is enhanced along the flow (x-) direction only (cf. 
Eqs. © and ©). This property is nicely reflected in the q- 
dependence of the growth rate, obtained from Eq. (fTTI) . In¬ 
deed, we And that S has no effect on the dependence of the real 
part of a upon qy, provided that qx = 0. In other words, the 
excitation and growth of all the modes of the type q = (0, g^) 
is independent of the differential flow parameter. As shown 
in Fig. [St, a necessary condition for the excitation and growth 
(Re [a] > 0) of these longitudinal modes is that e > r/Bf, 
since all modes with e < qB^ are damped. In order to con¬ 
trol mode selection, we first tune e in such a way as to obtain 
a positive growth rate for a longitudinal mode (see Fig. |S>), 
we then modify the growth rate associated with the transverse 
modes via a variation of S (see Fig. ©). As shown in Fig. 
a, for S = 0.7, the maximum growth rate for both modes 
becomes identical. For this set of parameters, both modes co¬ 
exist (see Fig. [tt,d). 

We numerically compute the stability diagram using the char¬ 
acteristic dispersion relation in Eq. (fT3l) . In Figure © we 
show a stability diagram of the modes, computed at param¬ 
eters 5 = 0 and e = 3.6. Modes within the grey area are 
unstable (give rise to patterns) and can coexist or interact, 
while those outside are damped. To access the stability dia¬ 
gram in (e, 6) plane, we consider modes with wavenumbers 
that fit into the typical size of our simulation domain {qx = 
27Tnx/Lx, qy = 2'Kny/Ly) and compute regions in (e, S) pa¬ 
rameter space where the growth rate Re [a] > 0. This is 
shown in Figur©?. In the regions designated by the ratio of 
the wavenumbers, the corresponding modes can coexist and 
interact. For instance the modes 4:2 corresponds to the ob¬ 
served pattern in Fig.©, while the mode 4:1 corresponds to 
the one in Fig.©. All other mode interaction can be obtained 
by choosing the corresponding parameters S and e accord¬ 
ingly. We expect that a wide variety of complex dynamics 
and spatial resonances can be triggered by the interaction of 
these steady state modes. 

In summary, the effect of shear flow on mode selection and 
length scale of patterns in nonlinear media is investigated. 
The mode selection under shear flow depends on the differ- 
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ential flow parameter of the participating species. Transverse 
and longitudinal modes are found to be selected depending on 
the differences in the flow rates of the participating species. 
Moreover, the effect of a heterogeneous shear on the underly¬ 
ing length scale of the patterns is accounted for via a simple 
scaling ansatz, which becomes exact at homogeneous shear. 
The predictions of the model regarding mode selection, length 
scale and the stability diagram are all in line with numerical 
simulation results. 


FIG. 4: (Color online) Plot showing dispersion relations at different 
differential flow parameters J (a) mode (b) comparison of the 
and qy mode at <5 = 0.7. 
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FIG. 5; (Color online) Stability diagram of excited modes, (a) In the 
qx — gyplane for 5 = 0.0 and e = 3.6. Beyond the grey domain, 
all the modes are damped, (b) In the (e, 5) parameter space showing 
the region for the wavenumbers (n^ : Uy) oi the modes. The system 
parameters are chosen as p = 2.1988, r = 2.744 and 7 = 0.001. 
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